function generate(v0,vth)
implicit none
real*8 :: generate, v0, vth
real*8, parameter :: pi= 3.14159265358979
real*8 :: sotto, vmax, fmax,v,ff,r,erf,ranf

sotto=v0*vth*sqrt(pi)/2.d0*(1-erf(-v0/vth))+vth**2/2.d0* exp(-v0**2/vth**2)
vmax=.5d0*(v0+sqrt(2.d0*vth**2+v0**2))
fmax=vmax*exp(-(vmax-v0)**2/vth**2)/sotto

1 v=(2.d0*vth+max(v0,0.d0))*ranf()
r=ranf()
ff=v*exp(-(v-v0)**2/vth**2)/sotto/fmax;
if(r<=ff) then 
 generate=v
else
 goto 1
end if
end function generate

